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Abstract 

, An extension of the H-iheorem for dissipative particle dynamics (DPD) to the case of a multi- 

component fluid is made. Detailed balance and an additional _ff-theorem are proved for an energy- 
conserving version of the DPD algorithm. The implications of these results for the statistical me- 
chanics of the method are discussed. 
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I. INTRODUCTION 



Interest in the rheological and dynamical properties of complex fluids [Q over the past decade has seen the intro- 
duction of several new techniques for their simulation on mesoscopic length scales. These methods include lattice gas 
automata (LGA), lattice Boltzmann equation (LBE) and dissipative particle dynamics (DPD) |^|]. The aim of this 
C letter is to explore some of the statistical mechanical properties of the rapidly evolving DPD model, and to extend 
these results in order to keep pace with new developments being made in the algorithms. 

The DPD method was originally introduced by Hoogerbrugge and Koelman || as a discrete time algorithm; this 
was subsequently modified and reinterpreted as a discrete time approximation to an underlying system obeying 
Langevin dynamics (with momentum conservation) by Espanol and Warren Q, in order to guarantee the existence 
of a Gibbsian (specifically a canonical) equilibrium state. Applications of the model include colloidal suspensions ||, 
polymer suspensions || and binary mixtures A dynamical theory has been presented for the continuous 

time limit and the equilibrium for finite time-step has been investigated [p"l| . 

We briefly describe the implementation of the DPD method. The system consists of a set of N discrete particles 
which move in continuous space and in discrete time-steps, the interval between which may be reduced to being 
infinitesimal. At each time-step St, the particles' momenta are updated by a momentum-conserving interaction 
with each particle inside a neighbourhood of radius Rq. This interaction includes three distinct forces, which can 
be described as conservative F c , dissipative F D and random F R . Between each tick of the clock, the particles all 
propagate freely according to their velocities. In the limit of continuous time, the DPD equations of motion are most 
effectively described in terms of the following stochastic differential equations: 
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q | for each particle labelled by the subscript i. The forces take the following forms: 

" ' n = (3) 

J m arij 

F § = -J w r>(rij)[eij ■ v SJ ]ey (4) 

Fy = <nu R (ry)eyCij (5) 

where <p is a potential energy, = — Vj is the relative separation vector and is the unit vector in the direction 
of Tij; for simplicity, all particles are assumed here to have the same mass m. The functions w D (rij) and w R (rij) are 



*email:c.marshl@physics.ox.ac.uk 

t email : coveney @ Cambridge .scr.slb.com 



1 



weighting functions which limit the action of the dissipative and random forces to a finite range Rq. The random 
elements Qj are Gaussian white noise with zero mean: Qj — 0. They are uncorrelated for different pairs of particles 
and for different times: Oj(t)(ki(t') = (Sik.Sji + 5u6jk)5(t — t'). It should be noted that these forces conserve both 
momentum and angular momentum but not energy. 

After a brief description of the statistical mechanical concepts involved, we present an ii-theorem for the multi- 
component DPD fluid. We then briefly describe the energy conserving DPD model and examine detailed balance and 
an ii-theorem in this context. 



II. ON DETAILED BALANCE AND ^-THEOREMS 



It is important at the outset to explain the significance and importance of the statistical mechanical properties of 
detailed balance and the so-called ii-theorem. Detailed balance is known to be a sufficient, but not necessary, condition 
which ensures that a Gibbsian equilibrium state exists in the ensemble representation of a dynamical system. It is 
possible for systems not satisfying detailed balance to exhibit equilibrium states; however, characterising these is then 
a much harder task. The virtue of the modifications made by Espanol and Warren to the original DPD algorithm is 
that, in the limit of continuous time, the iV-body DPD system then satisfies the detailed balance condition, thereby 
guaranteeing the existence of a well defined equilibrium state. 

In the literature, at least two separate kinds of "ii-theorem" can be distinguished. First, any Markov chain or 
process which has an equilibrium distribution will have an ii-theorem associated with it, in the sense of possessing 
a Lyapounov function that changes monotonically with time. Indeed, a whole class of Lyapounov functions achieve 
this; the class is defined as the expectation of any convex function of the relative-to-equilibrium probability density. 
The proof of such ii-theorems follows directly from the linear equation (also referred to as the "master equation") 
for the iV-body distribution function. Detailed balance plays a role here, in that it specifies what the equilibrium 
distribution is; this information is needed to write down the Lyapounov function. 

However, the arguably more famous Boltzmann ii-theorem is quite a different notion, and is of much more restricted 
validity. Boltzmann's H function is defined in terms of the one-body distribution function, and its time-monotonicity 
can only be derived if we know the kinetic equation obeyed by this reduced distribution. Moreover, this equation is 
non-linear so that the choice of H is now much more restricted; for the Boltzmann equation itself, only the expectation 
of the logarithm of this reduced probability distribution suffices. In proving Boltzmann's ii-theorem, use is made of 
the property of detailed balance. 

Note in passing that time-symmetry is a stronger property than detailed balance; the former implies the latter, 
but is not implied by it. Thus, detailed balance is obeyed by both Newton's equations of motion and by dissipative 
particle dynamics, although the former is time-symmetric while the latter is not. 

The existence of an ii-theorem for a given system can be used to check on the numerical stability of any algorithm 
implemented to simulate it; numerical instabilities which lead to non-monotonicity of the H- functional concerned can 
then be precluded. The issue of the existence of detailed balance and related ii-theorems is thus clearly of importance 
for the various mesoscale modelling and simulation techniques. 

By contrast with (continuous time) DPD, virtually all interacting lattice-gas and lattice-Boltzmann models have 
no known detailed balance condition; therefore, their equilibrium states are generally unknown, while the lack of 
any associated ii-theorems makes the real-valued lattice-Boltzmann methods, in particular, susceptible to poorly 
understood numerical instabilities. Indeed, because detailed balance is not satisfied in such models, it makes their 
theoretical analysis by standard methods of non-equilibrium statistical mechanics well nigh impossible. 



III. if-THEOREM FOR MULTICOMPONENT, ISOTHERMAL DPD 



Detailed balance and an ii-theorem (of the first kind mentioned in the preceding section, i.e. for the full TV-body 
distribution) for the single component DPD fluid have already been demonstrated M. The proof of detailed balance 
for general DPD models of interacting multi-component fluids has also been derived Here, we aim to extend this 
form of ii-theorem to the case of a multi-component fluid, which includes the case of binary immiscible fluids 0,8|. 

It has been demonstrated (l3) that the evolution equation for the iV-particle distribution function is the Fokker- 
Planck equation for the multi-component fluid: dfP = C MC P, where the multi-component Fokker-Planck operator 
£ MC is defined as: 
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where a and (3 are sums over different types of particles and i a and jp sum over all particles of each type. The 
parameter 9 is defined as 9 = mo~ 2 /2j. The relevant H-functional for the multi-component case is a simple extension 
of that in the single component case. As expected for an isothermal system, it is just the expectation of the associated 
free energy < U — OS >, where U is the internal energy, 9 is the equilibrium temperature, S is the global entropy, 
and the expectation is taken using the full iV-particle distribution function, P: 
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Using the time evolution operator for the multi-component system (ra), it is possible to show that: 
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It is then apparent that the time derivative of the functional T is the sum of negative definite terms, and therefore 
that the functional itself is monotonically decreasing in time. The appropriate equilibrium distribution for the multi- 
component system occurs when this functional stops decreasing. It is easy to show that this occurs when it reaches 
the Gibbsian distribution for the associated conservative system, i.e. as if the dissipative and random forces were not 
present: 



• exp 
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(9) 



Z MC being the multi-component canonical partition function, defined in the normal way. 



IV. ENERGY CONSERVING DPD 

An energy conserving version of DPD has recently been presented by Espanol |L2| . This involves the introduction 
of an internal energy variable for each DPD particle (now interpreted as a cluster of atoms or molecules, into which 
the dissipated energy is assumed to flow). There is an entropy s(ej) which needs to be specified in order to describe a 
given system, and the temperature is defined in the usual thermodynamic way as 9% = (dsi/dei)^ 1 . It is then possible 
to formulate a new DPD algorithm which conserves the total energy of the system, as well as momentum and angular 
momentum [O. An appropriate set of stochastic differential equations is: 
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(12) 



Here the functions A D {rij ) and A R (r^ ) are additional weighting functions for what can be interpreted as the conduction 
and random heat flux terms respectively, while K%j and are their strengths; Q- are random elements which are 

uncorrelated to the elements Qj and have zero mean Qj = 0. They are uncorrelated for different times and different 

pairs of particles and are anti-symmetric: CyWCjwC^) = ($ikO~ji ~ Su5jk)5(t — t'). We also note that the strengths of 
the random and dissipative forces {a and 7) can now, in general, vary for different particle pairs. 
Following the original derivation |T^|, we make the additional assumptions: 



Al(r) = A R (r), 



w 2 R (r) = w D (r). 



(13) 
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which mean that the Fokker-Planck equation for the evolution of the TV-particle distribution function P, can be 
written as: 



d t P = [jC c + C VH + C HC ] P 
where the operators on the right hand side are defined as follows: 
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(17) 
(18) 



where the subscripts C, VH and HC refer to the Conservative, Viscous Heating and -ffeai Conduction terms 
respectively. 



V. DETAILED BALANCE FOR ENERGY-CONSERVING DPD 

If the evolution operator for a system is C and we designate its adjoint by the operator C\ then the detailed balance 
constraint p4[ can be written in the following way: 



£*Peqm^P — Peqm^ '■P 



(19) 



where the superscript e indicates that all variables that are odd under time reversal are to have their signs reversed. In 
the case of DPD, this means that the velocities attract an additional minus sign. The function ip can be any function 
of the phase space variables. 

The appropriate operators for the energy-conserving version of DPD are: 
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We can then show that: 
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and 
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where we choose the strengths of the dissipative and random forces to satisfy the following relations: 
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It is therefore apparent that: 

[£c + £VH + jC-Hc]Peqm(p = Peqm[£>G + ^VH + £hc]*P 



(25) 
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and therefore that the energy-conserving DPD algorithm satisfies detailed balance. The equilibrium distribution 
associated with this detailed balance condition satisfies the following relations: 
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Therefore the equilibrium distribution consistent with [Cc + C-vh + C-Hc\Peqm — has the following form: 
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where Z EC is the normalization constant. 



VI. ^-THEOREM FOR ENERGY-CONSERVING DPD 

The ff-theorem for energy-conserving DPD can now be formulated. We first define the following iJ-functional: 



S[P(T)](t) = J dTP(T,t) 



J> -lnP(r,i) 



(31) 



Then, with the aid of the appropriate Fokker-Planck equation (|1J), it is possible to show that the time-evolution of 
this functional is: 
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We note that the time derivative of the functional S consists of sums of two types of terms, each of which is positive 
definite. Therefore S^PfT)] is monotonically increasing in time and the equilibrium is reached when this time evolution 
stops. It is easy to show that this can only be achieved when the equilibrium distribution satisfies the following 
relations: 
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and therefore that the equilibrium distribution consistent with [Cc + £vh + C-Hc\Peqm — will be: 



P, 



eqm 



■ exp 



^2 Si (Ei) 



(35) 



This is naturally the same equilibrium distribution already recognised as being a stationary point of the Fokker-Planck 
evolution eqns ( p8| ) and (^9|). The TV-body P-theorem provides additional information because it guarantees that the 
system approaches this equilibrium state monotonically. 

The ^-functional (|l|) may be interpreted as the total entropy of the system. We see that the first term in the 
functional represents the microscopic entropy of each DPD particle while the second term represents the normal 
macroscopic entropy —Pin P. That the relevant functional is in this case the total system entropy, rather than a free 
energy, could be expected from the fact that the system is now energy conserving. 

For notational ease, the detailed balance property and P-theorem for DPD have been presented for the single 
component case. However, their extension to the multi-component case can be achieved in a similar manner to the 
isothermal P-theorem result presented here. 



VII. CONCLUSIONS 



We have shown that the desirable statistical mechanical properties of detailed balance and the existence of P- 
theorems may be extended to general multi-component interacting DPD systems, whether maintained at constant 
temperature or at constant energy. Of course, such properties are rigorously valid in the continuous time limit, and 
are only approximately true for discrete-time implementations of these algorithms. The approximations improve as 
the size of the time-step is decreased. Detailed balance makes possible the theoretical analysis of models based on 
the DPD equations of motion, while the P-theorem provides a means to control numerical instabilities in computer 
simulations. 
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